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Abstract 



^vq We have implemented open boundary conditions into the ANTARES code to increase the realism in our simulations of 
stellar surface convection. Even though we greatly benefit from the high accuracy due to our 5th-order numerical scheme 
(-^ (WEN05), the broader stencils needed for the numerical scheme complicate the implementation of boundary conditions. 
Moreover, one-dimensional starting models may not provide enough data to set up the simulation domain large enough 
in each direction, or the initial data may lead to numerical instabilities. We show how to overcome these difficulties 

CO and demonstrate to what extent numerical simulations of stellar surface convection are sensitive to the numerical setup 
and the boundary conditions. A sloppy choice of parameters contained in the boundary conditions can have a severe 

p/| impact. Only if the initial model, the extent and position of the simulation box, and the parameters from the boundary 
conditions are chosen adequately, numerical simulations of stellar surface convection are (physically) meaningful and 
• realistic throughout the entire simulation domain. 

Q^JCeywords: hydrodynamics - methods: numerical - stars - Sun: granulation - convection 
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1. Introduction 

Nowadays, realistic simulations of stellar surface con- 
vection are a mature tool of computational astrophysics, 
and find a wide field of applicat ions. An exten s ive re - 
vi ew on this subject is given in INordlund etaD fl2009h . 
In lBeeck et al.l (|2012l ). the usefulness and reliability of re- 
sults obtained from simulations of solar surface convection 
With the codes C05BOLD, MURaM and Stagger is shown. 
Even though numerical methods and the setup of the sim- 
ulations in this comparison are quite different, the overall 
stratification and other basic properties of the numerical 
model are similar. 

For any numerical simulation which requires the solu- 
tion of partial differential equations, the choice of the nu- 
merical method, the simulation domain, and appropriate 
boundary conditions are of major importance. For the 
simulation of convection at the surface of solar-type stars, 
the top boundary of the simulation domain is in the up- 
per photosphere, whereas the bottom boundary is situated 
deeply inside the convective envelope. However, due to the 
mixed parabolic-hyperbolic nature of the governing equa- 
tions the boundary conditions influence the solution in the 
whole simulation domain. 

The ANTARES code has originally been devel- 
oped at the University of Vienna for the simula- 



tion of surface c o nvecti on (jMuthsam et al.l . 120071 [2010a; 
Lemmerer et all 120131 ). Recently, the code was 
also applied to many other astrophysical problems 



(Kupka ct al., 2009; Muthsam et al., 2010b; Kupka et al 
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2012; Happenho fer et al.. 12013c iZaussinger and Spruitl . 
2013J; iMundprecht et all 20131) . In some of these appli- 
cations pro blems arose due to t he use of closed boundar y 
conditions (|Kupka et al.l , 120091: IMundprecht et all 120131 ). 
Shock fronts were reflected at the top boundary, leading 
to instabilities in the simulation results. Boundary condi- 
tions which transmit shock waves rather than reflect them 
are needed to continue these simulations over a sufficiently 
long time interval. 

iRobinson et al.l J2003h . iKupka and Robinson! (j2007l) 
and iKupkal (|2009allbh discussed, amongst others, the ef- 
fect of boundary conditions on statistical properties of the 
flow. Unrealistic boundary conditions can lead to unphys- 
ical flow patterns in a huge part of the simulation box. 
When the bottom boundary is situated in the convective 
region of the star, which usually is the case, e.g., for so- 
lar surface convection simulations, most of the energy is 
transported in the lower part of the simulation box by ad- 
vection of enthalpy and kinetic motions. Therefore, free 
in- and outflow into the simulation domain is necessary 
to avoid an artificial flow pattern and artificial heating in 
the lower part of the simulation domain. The entropy of 
the inflowing material is unknown and has to be specified 
somehow. At the top boundary, shock fronts should not 
be reflected, but be transmitted. 
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Even though a lot of literature exists on stellar in Section [5j 



surface convection simulations ( e.g.. iNordlundl . [1982 



Frevtag et all l2012t IVogler et all [2005), a dedicated in- 
vestigation on the effect of the setup of boundary condi- 
tions on the simulation results appears to be missing. This 
paper gives a precise description of the numerical setup 
of our simulations with ANTARES and demonstrates the 
sensitivity of the simulations to changes in the choice of 
parameters which are introduced by specifying the bound- 
ary conditions. 

As we can argu e when considering the study of 
Asplund et all (|2000h . meaningful tests of boundary con- 
ditions can ultimately only be done in three dimensions, 
with sufficient resolution and long relaxation time. The 
two-dimensional geometry leads to systematic differences 
in the mean stratification such that, for instance, a dif- 
ferent inflow entropy is required to obtain the same total 
flux at the stellar surface. This makes this investigation 
very expensive in terms of computation time as well as in 
"real" time. 

The results of the present paper provide a nece ssary 



extension to comparisons shown in iKupkal (|2009ah and 
( 20121 ). where simulation data from very dif- 



Beeck etal 



ferent simulation runs were compared. Systematically 
testing every parameter combination in the setup of a 
numerical simulation, however, is simply not affordable, 
whence we have to rely on tests of the overall bound- 
ary condition with only partial coverage of the parameter 
space available. In this paper, we show results from several 
working and non-working combinations and estimate the 
impact of inappropriately designed boundary conditions. 

The remainder of the paper is organised as follows. In 
Section 12.11 we show the superiority of the WENO algo- 
rithm used in ANTARES compared to ordinary first or 
second order schemes. The gain in accuracy justifies the 
broader stencils needed for any high-order method. On 
the other hand, the design of boundary conditions gets 
more complicated since they extend over several vertical 
layers. We have implemented several boundary conditions 
also used in the aforementioned codes paying special at- 
tention to the requirements of our high-order numerical 
scheme. Details are given in Section 12.21 The numerical 
setup of our simulations, i.e. the starting model and the 
simulation box size, is described thoroughly in Section [23] 

We then investigate the dependence of the boundary 
conditions on parameters like the amount of energy flow- 
ing into the domain at the bottom boundary and com- 
pare them, taking statistical properties of the flow into 
account, in Sections 13.11 and 13.21 We show the difficul- 
ties arising in the setup of the simulation, its dependence 
on the initial model (Section 13. 3p and on the grid reso- 
lution (Section 13.4)) . and emphasize the importance of a 
careful choice of parameters. In Section 13.51 we describe 
differences between simulations with open and with closed 
boundary conditions. Finally, we compare two- and three- 
dimensional models in Section 13.61 A discussion of the re- 
sults can be found in Section @] followed by our conclusions 



2. Simulation Setup 

In this section, we describe the implementation of sev- 
eral boundary conditions in the code ANTARES as well 
as other modifications to the code we have performed in 
this context. 

In this document, the ANTARES convention for spatial 
coordinates is used, i.e. the x direction points inwards into 
the star and u is the vertical component of the velocity 
vector (u, v, w) T . The grid is assumed to be equidistant in 
every direction. 

2.1. Some Comments on the Numerical Method 

For the spatial discretisation of the advective part of the 
Navier-Stokes equations, ANTARES uses the weighted es- 
sentially nom^o^dlla^ory (WENO) a lgori t hm as described, 
e.g., in IJiang and Shu] (Il996h and [Sh] ([20031) . WENO 
schemes exist for any order of accuracy. In ANTARES, 
we use the fifth or der variant which will b e called WEN05 
in the following ( Muthsam et all l2010al) . For the tests 
discussed in this subsection , flux splitting as described in 
Donat and Marquina ( 1996[) is not necessary, even though 
it is implemented in ANTARES and used for the simula- 
tions presented in the later sections of this paper. 

The advantages of the WEN05 method are its high or- 
der of accuracy and its shock-capturing ability. This al- 
lows accurate modelling of surface convection where strong 
shock fronts are ubiquitous in the downflo ws, whereas the 
granu les themselves are rather smooth (cf. iNordlund et al. 
2009). Even when used in combination with lower order 



methods, e.g. for the time integration, which theoretically 
degrades the overall order of the method, the WEN05 
scheme is superior to classical schemes or other ENO- 
type schemes of lower order, especially when the latter 



are combine d with artificial diffusivities (jMuthsam et al 
2007ll2010ah . 



We demonstrate the efficiency of the WEN05 method 
by solving the advection equation 



dt dx 



(1) 



for t > and x € [0, 1] with periodic boundary conditions. 
With the initial condition 



(j){x,0) = 1 + 0.1sin(27ra;) 



(2) 



the analytical solution stays smooth for all times. There- 
fore, this is an appropriate test case for determining the 
empirical order of accuracy and the error constants of a 
method. In these test calculations, the Courant number is 
fixed to 0.1. 

From Figure [1] we deduce that the empirical order of 
accuracy of the WENO 5 algorithm toge ther with a sec- 
ond o rder time integration such as TVD2 ([Shu and Osher , 
1988ft is two which is also obtained with the Lax-Friedrichs 
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Figure 1: Empirical order of accuracy of several schemes when solv- 
ing iTiJ and f2t . The thin lines indicating first and second order 
convergence are given for guidance. 



Figure 2: Empirical order of accuracy of several schemes when solv- 
ing Q and JS}. The thin lines indicating square root and first order 
convergence are given for guidance. 



method (cf. IStrikwerda 19891) . Nevertheless, the error 
is two magnitudes smaller than with the Lax-Friedrichs 
method except for the small region shown to the right 
where the resolution is very low. There, we observe third 
to fourth order convergence of the WEN05 method since 
in this region, the spatial error dominates over the tempo- 
ral error and the higher-order convergence of the WEN05 
scheme is visible. The error of any first order method as, 
e.g., the upwind method (cf. again IStrikwerda 1989), is 
larger by several magnitudes. A very high amount of grid 
points is required for them to reach an acceptable error 
size. 

Typically, the magnitude of the (relative) mean error 
in our simulations is of the order 10 -3 . With the Lax- 
Friedrichs method, we would need about six times as many 
grid points compared to the WEN05 method to obtain the 
same accuracy. This means that in three dimensions the 
number of grid points must be more than 200 times higher 
than with the WEN05 method. Furthermore, this accu- 
racy is only achieved with a much smaller time step size. 
In the end, for smooth flows we have to expect that the 
right-hand side has to be evaluated over a 1000 times more 
often when using the Lax-Friedrichs method in compari- 
son with the combination of TVD2 with WEN05 even at 
the moderate resolution typical for astrophysical simula- 
tions (~ 4 grid points per smallest feature) to achieve the 
same accuracy. This can be summarised by the combined 
scheme of WEN05 with TVD2 having a much smaller er- 
ror constant than the Lax-Friedrichs scheme. 

Therefore, the WEN05 method is not only superior in 
terms of shock-handling, but also in terms of efficiency. In 
contrast, the disadvantages — higher computational costs 
per grid point and broader stencils — are negligible. 

Given discontinuous initial data, 



the convergence order is restricted by the smoothness of 
the (weak) solution. In Figure [5J the order of all three 
methods is worse than linear. Nevertheless, WEN05 
with TVD2 is again much more efficient than the Lax- 
Friedrichs method: for an error of magnitude 0.1, about 
six times as many grid points are needed by the Lax- 
Friedrichs method compared to WEN05. Again, the first- 
order accurate upwind method performs worst and shows 
the largest error except for very coarse resolution. 

In conclusion, for both smooth and discontinuous solu- 
tions the WEN05 method combined with TVD2 time inte- 
gration is more efficient and more accurate than standard 
low order methods. This justifies the additional efforts re- 
quired for its implementation such as adapting boundary 
conditions to its wider stencils. 



Fro m the comparisons done in iMuthsam et all (|2007f) 



and in IMuthsam et all (|2010al) . it follows that WEN05 
without artificial diffusivities yields also much more accu- 
rate results than other high-order methods which, how- 
ever, always require such stabilisations. Also in sim- 
ulations based on the WEN05 numerical scheme unre- 
solved sc ales can be modelle d with a subgrid scale mo- 
del (e.g., Smaeorinskv . 19631) . Anyhow this was not nec- 
essary for the simulations presented in this paper, since 
their numerical resolution was rather high and no partic- 
ular stability problems have been encountere d (implicit 
LES or LLES approach, iGrinstein et all 120071 ). We also 



note that the appropriate treatment of diffusive and vis- 
cous terms in the context of WEN05 is discussed in detail 



Happenhofer et all (|2013l) . 



0(3,0) 



1, if 0.1 < X < 0.3, 
0, else, 



(3) 



When applying the WENO method to systems of con- 
servation laws, the state variables must be transformed 
into the eigenstate which increases the computation time. 
On the other hand, methods where no transformation is 
needed are less accurate and artificial diffusivities are nec- 
essary to stabilise the solution, e.g. arou nd shock fronts , 
which at the bottom line is less efficient (Muthsam et al 
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2007 L l2010ah . 



2.2. Boundary Conditions for Stellar Convection Simula- 
tions 

In simulations of surface granulation of cool stars, the 
top boundary is typically placed above the photosphere 
whereas the b ottom boundary is deep inside the convec- 
tive envelope dstein and Nordlundl, 1998 : Robinson et al 



20031: IWedemever Bohm et all 120041: iFreytag et all 12012 
Voder et all l2005HMuthsam et alll2007ll2010;i 7Numer- 



ical boundary conditions either force the fluid to stop its 
vertical motion and move horizontally without friction (the 
so called "slip" boundary condition) and hence not cross 
the vertical boundary, or they allow (constant) in- and 
outflow. In the case of an inflow, the density and the en- 
ergy of the inflowing material must be specified. 

There is no physical reason why the vertical convec- 
tive motions should stop at the boundary of the simu- 
lation domain as it is forced by closed boundary condi- 
tions. To increase the physical realism of the simulation it 
is mandatory to allow free in- and outflow at the bound- 
aries. On the other hand, implementation and stability of 
open boundary conditions are non-trivial, as will be shown 
in the following. In particular, the number of boundary 
layers must correspond to the width of the stencils of the 
numerical method used. 

For solar-like main-sequence stars, the characteristic 
scale of surface convection is small compared to the stel- 
lar radius and, hence, the "box-in-a-star" approximation 
is valid for the simulation of surface convection. If the 
simulation box is broad enough, periodic boundary con- 
ditions can be used on the horizontal boundaries of the 
domain without any n egative impact on the simulation 
( Robinson et all 120031) . 



2.2.1. Closed Boundaries 

Slip (Stress-free) boundary conditions are characterised 
by setting the vertical velocity component and the vertical 
derivative of the horizontal velocities to 0, i.e. 



„ dv dw 
dx ox 



(4a) 



Furthermore, since the bottom boundary can transport 
energy only by radiation, an artificial source term must be 
introduced to feed the required amount of energy into the 
simulation domain. This is done by modifying the thermal 
conductivity k in the lowermost layers of the simulation 
domain such that 



8T - P - 7* 
^modified q — — ^^eff' 



(4b) 



a nd imposing ^ = 0. Th e det ailed setup is described, e.g. 



Robinson et ail (j2003l ) and iMuthsam et al.l (|2010al ). 



Closed boundary conditions have the advantage of easy 
implementation and high stability, at least for subsonic 



flows. On the other hand, they reflect shocks and dis- 
turb the velocity field in an undesirable way. In the con- 
text of simulations of stellar surface convection, their un- 
favourable effect concerning reflection of waves and sta- 
tistical pr operties of the flow is discussed in detail, f or in- 



stance, in 



Robinson et al.1 (|2003luKupka et al.l (|2009l) . and 



Kupkal ( 2009a ). Nevertheless, there are applications where 
the use of closed boundary conditions is not problematic, 
such as the nume rical simulation of the superadiabatic 
layer of the Sun dKim and Chan . 1998 : Robinson et al 



20031: iKupka and Robinsonl . l2007t iKupkal . l2009a|) 



2.2.2. Open Top Boundary Conditions 

T he boundary c o nditio ns presented in Cheund ( 2006j) 
and Trampedach] (1997 ) and also those discussed in 
Frevtag et al.l (|2012l) are not directly applicable to AN- 
TARES due to the broader stencil of the WEN05 method 
and the equidistant grid. Therefore, we had to adapt their 
approach to our setup. 

For the top boundary layers, we assume a hydrostatic 
and isoenergetic, i.e., in practice, isothermal stratification. 
As a difference to previous implementations, in the hydro- 
static equilibrium we want to include the turbulent pres- 
sure since it is quite large in the photosphere. The veloc- 
ities should be constantly extrapolated. Generally speak- 
ing, the boundary conditions should be designed in such 
a way that the numerical setup does not influence their 
behaviour. 

The vertical and horizontal velocities u, v, and w in the 
boundary layers are set according to 



du dv dw 
dx dx dx 



(5a) 



For a high-order method, more than one boundary layer 
is needed. Therefore, we add as many additional nodes on 
the top of the first boundary layer as necessary for our 
scheme. In the case of WEN05, three layers are required. 

Let be the index of the innermost boundary layer and 1 
the one of the outermost domain layer. The mean density 
in the innermost boundary layer (p)o is advanced in time 
by an averaged midpoint rule centred at the outermost cell 
interface of the domain, i.e. 

(p)(«+ 1 ) = (p)W-^( pu) (»), ( 5b ) 

L\X 2 

where T stg is the time step size of the current Runge- 
Kutta stage and Aa; is the (constant) vertical grid spacing. 
(pu) " , the mean density flux at the interface between the 

2 

layers and 1, can be computed by the spatial integration 
scheme of its dynamical equation given enough boundary 
nodes. Now, the density profile from the outermost do- 
main layer is just copied to the innermost boundary layer 

scaled by the ratio 



(p) ( r r +1) 



p^ 1 HQ,V,z) = y^ffi* 1 \l, V ,z). (5c) 
(P)i 



4 



The specific internal energy (e)o in all boundary layers 
is assumed to be constant. It is initialised by the mean 
specific internal energy in the outermost domain layer (e) 1 
and is furthermore advanced in time by 



(4 n+1) 



, -(1.0-6){e)™+S(e)r X1 - (5d) 

ICheund (2006) set d constantly to 10~ 3 , but we reset (e) 
only at the end of each Runge-Kutta step with 5 given by 
the formula 



5 = 



r • (Vsndh 



(5e) 



where r is the time step of one Runge-Kutta integration, 
(^snd)i is the mean sound speed in layer 1 and (H)\ the 
mean pressure scale height in that layer. 5 £ [0, 1] controls 
the time scale on which the heat content of the boundary 
layers varies. With the parameter Cf, this time scale can 
be changed: higher values of Cf lead to slower changes, 
whereas with small values of Cf the boundary layers relax 
much faster to the stratification of the inner cells. Since 
Cf has to be adjusted anyway, it is sufficient to use the gas 
pressure to compute (H)\ for the approximate time scale 
underlying (|5e|) . 

The idea behind equation (|5e|) is to make the procedure 
independent of different time integration methods and spa- 
tial resolution. With a constant value of S, the time scale of 
changing (e)o would be different, if a two- or a three-stage 
Runge-Kutta method is used, or if the Courant number 
is changed. Furthermore, we choose <5 independent of the 
grid spacing such that the time scale does not change, if 
the grid is modified. Modifying Cf allows to control the 
stability and the stratification of the boundary layers. 

The density in the layers above index is set according 

to 



dip + Pt» 
dx 



-P9, (5f) 



where p = p rad + p gas and p tu rb = p(u 
bulent pressure. By the chain rule, 



{it}) is the tur- 



0(p + p tv 



<9(P+Pturb) dp 



dx 



dp 



dp , , ,, 2 \ dp 
dx-\dp- + {u ~ {u)) )d-x 

(5g) 

For an ideal gas, constant e implies an isothermal strat- 
ification. The slight fluctuations of T, evaluated from the 
realistic equation of state, therefore show the deviation 
from the ideal gas equation. Numerical tests of this con- 
dition show that it is reasonable to assume the simplifica- 
tion of an ideal gas in the boundary layers, considering the 
other simplifications like constant velocity and constant e. 
In this case, = (7 — 1) e is constant in the boundary 
layer. 

Since (u — (u)) 2 in the boundary layer is constant in 
the vertical direction due to (|5"aj) . the value c p := + 

(u — (u)) 2 is constant for each horizontal grid index (y, z). 



Therefore, we can integrate equation fl5g| ) analytically and 
obtain 



p{hV,z) = p(0,y,z) exp 



— \i\ Ax ■ g 



(5h) 



where i indexes all boundary layers above layer 0. With 
(e)o and the velocities calculated by (|5aj) , all values in all 
boundary layers can be calculated. 

2.2.3. Open Bottom Boundary Conditions 

At the bottom boundary, we set the velocities according 

to 



du 
dx 



dv 
dx 



dw 
dx 



= 



(G) 



to allow constant in- and outflow. In the literature, the 
density (of an inflow) usually is set to keep the total 
mass of the model unchanged. Therefore, only an ad- 
ditional condition for the energy is needed. Since the 
bottom boundary in surface convection simulations is sit- 
uated in the adiabatic region, a constant entropy value 
ca n be assumed for infl o ws. T his idea was first presented 
in iNordlund and Steinl dl990l) and is w idely used today 



Frevtag et al.U2012t IVogler et all 120051) . 



There are two approaches how the amount of energy of 
a n inflow at the bottom bou ndar y is specified technica lly: 
in lNordlund and Steinl (|l990h and lFrevtag et al.l ([20121 ) for 



instance, an e ntrop y value is prescribed, whereas, e.g., in 
Vogler et al. (|2005l) . the specific internal energy is pre- 
scribed. Since the correct values are not known from 
within the problem in both cases, they must be found 
iteratively or be taken from some external source. For 
small changes, the relation between e and S is nearly lin- 
ear, so that it is not important to which of these variables 
the boundary conditions arc applied to. Nevertheless, e in- 
creases with depth, whereas S in the deep layers of the con- 
vection zone should approach an (almost) constant value. 
Therefore, prescribing an entropy value has the advantage 
of being essentially depth-independent, whereas values for 
the inflowing energy must vary with box depth. 

The density in the boundary layers in both cases is set 
such that the total mass contained in the simulation do- 
main stays unchanged. 

Since in a relaxed simulation, the radiative flux at the 
top of the domain corresponds to th e amount of energy 



top 01 tne domain corresponds to tne am ount ot energy 
flowing through the bottom boundary. IVogler et al. (2005) 
decided to directly correlate the two values on the Kelvin- 
Hclmholtz time scale. When the radiative flux at the top 
is too high, the amount of energy flowing in at the bottom 
is lowered which in principle leads to a lower radiative flux. 
This motivated the formula 



(n+l) 
inflow 



» 

"inflow 



1.0 



1.0 



tkh 



ptop 
J rad 



(7) 



5 



where = <jT^ b is the energy flux of the star accord- 
ing to the Stefan-Boltzmann law and r the time step size 
of one Runge-Kutta integration. The Kelvin-Helmholtz 
time scale tkh can be calculated by 



tkh 



J top Fad dydz 



(8) 



(cf. formula (19) a nd (20) inlVogler et all (120051) 1. 

On the contrary, iNordlund and Steinl (jl990T) decided to 
choose a fixed value for the entropy of the inflow at the bot- 
tom boundary, although it is not di scussed how in genera l 
this value of Si n fl ow is obtained. In iFrevtag et al.l (|2012l ). 
the mass and the energy of the inflowing material are cor- 
rected in two steps to keep the entropy of the inflow close 
to Sinflow and to reduce pressure fluctuations (see their 
equations (34) to (37)). They also do not discuss how the 
value Sinflow is actually obtained. 

We implemented both approaches which define either 
Sinflow or einflow into our code and will show results in 
Section 13.21 In those cases where we specify S; n flow at the 
bottom, we use an iterative correction procedure similar 
to ([7]). which is used for einflow- In the first variant, we 
only change the time scale and define 



o(n+l) _ 
inflow inflow 



1.0 



I.O- 



ts 



ptop 
r rad 



(9) 



where r$ is the time scale of the correction. In general, 
T s 7^ tkh since the correction procedure is applied to S and 
not to e as in (|7J| . In the limits r§ — > oo and F £f — > J* 1 *, we 
arrive at a constant value of Si n fl ow w hich is effectively the 
bound ary condition that was us ed bv lNordlund and Steinl 
(jl990l) and lFrevtag et al.l (j2012t) in the sense that S; n flow = 
constant. 

In the second variant, it is guaranteed that the total flux 
at the bottom boundary is close to F+. Here, 



inflow 



7 (") 
"mflo 



1.0+1- 



1.0- 



F 



bo I 



F* 



(10) 



where is the sum of the radiative, the kinetic, and 

the convective fluxes at the bottom of the domain. 

None of these correction mechanisms should be applied 
in the first minutes of the time evolution of a model, since 
both the radiative flux at the top and the total flux at 
the bottom require some time to reach the value which 
reflects the current numerical setup and thermal stratifi- 
cation. Therefore, we do not change S; n flow during the first 
five sound crossing times of a new model started from a 
one-dimensional stratification. 

Specifically, in the first correction step the density p and 
specific internal energy e of an inflow are changed such that 



S(pW e«) = Si 



inflow j 



(ii) 



where Sinflow is determined either with ([9]) or (fT0|) . En- 
tropy, internal energy, and dens ity of an outfl o w are not 
changed. On the other hand, in IFrevtag et al.1 (|201 2| ) the 
density and specific internal energy are reset by 



(i) , t yr(r 3 -i) 

P = P + CSchangeT ~p (.^inflow ~ O j 



tch&r pTl 

r 3 -i 



£ (1) = e + CSchangeT T [ 1 - 

^char 



(12a) 

(Sinflow — S) , 

(12b) 



(cf. (34) and (35) from IFrevtag et alj|2012» with the pa- 
rameter cschange controlling the time scale of the correction 
step (further quantities are explained below). In our im- 
plementation, however, the entropy of the inflowing mate- 
rial is fixed to S;^^* determined by © or (fill)) , therefore 
deleting the parameter cg c hangc from our formulation of the 
boundary conditions. In the end, keeping the gas pressure 
unchanged, we compute Sinflow := Sj"g^^ from either ([9} 
or (TTOl) and obtain and p^ from inverting (fTTj) using 
the equation of state. 

We then proceed with two ad ditional correction steps 
proposed bv IFrevtag et al.1 ( 2012 ) which avoid the genera- 
tion of unwanted large pressure fluctuations by the inflow 
at the open bottom boundary, and fix the mean mass flux 
over the boundary to 0. Thus, as in the second correction 
step from IFrevtag et al.l ( 2012 ) , density and specific inter- 
nal energy throughout the entire boundary layers are reset 
by 



p(2) = p(l) _|_ Cpcnangc 



— «p) - p) , (13a) 

tchar V snd 



.(1) 



CpchangcT ^— ( (p) ~ P) (13b) 

tchar 1 lP 

(cf. (36) and (37) from IFrevtag et aL I l2012l) . Here, r is the 
time step size of the (Runge-Kutta) integration and £ cna r 
is given by 



^char 



Ax 



(fsnd + |«|) 



(13c) 



(cf. (33) from lFrevtag et al.ll2012l ). The velocity of sound 
u sn d as well as the adiabatic coefficients Ti and T 3 are 
taken from the equation of state as a function of the uncor- 
rected values of p and e. The parameter cp cn ange governs 
the time scale on which pressure fluctuations are reduced. 
Finally, p and the vertical velocity u are modified by 



p(0 + (,)(«» _(,»), 



(p^u) 



(13d) 
(13e) 



G 



(cf. (38) and (39) from iFrevtag et alJlioil) to keep the 
total mass at the bottom boundary unaltered. Again, this 
correction is applied throughout the entire boundary lay- 
ers. 

As for the top boundary at the bottom three bound- 
ary layers are needed for the WENO scheme used in AN- 
TARES. After the procedure described above is applied 
to the innermost boundary layer at the bottom, the two 
underlying layers are filled by exponentially extrapolat- 
ing the density and by linearly extrapolating the specific 
internal energy (assuming adiabatic stratification instead 
would require to proceed as in Sect. 12.3.21 to recover the 
power law growth with depth, possibly with a correction 
for turbulent pressure as in f5j]) and fl5gl ), but we expect 
the difference to our simplified procedure along the two 
outer grid points to be acceptably small). 

When i mplementing the bou ndary conditions which are 
similar to Frevtag et al. ( 20121 ) , care must be taken in the 
formulation of mass conservation. Setting the mean of pu 
as a variable to does not enforce mass conservation in 
flux-based codes such as ANTARES. Instead, the corre- 
sponding cell boundary flux must be modified. The system 
is not overdetermined with this condition, since it is only 
another numerical enforcement of the (analytical) mass 
conservation requirement. 

Furthermore, as another available alternative to the 
open bottom boundary conditions specified by ([B]), © 
or (fT0]l . (jlll) and (|13p . we i mplemented t he fu ll version 
of boundary conditions from Vogler et al. (2005) into our 
code. For the present paper these are used only for the 
case of Model 2 in Section 13.21 In these boundary condi- 
tions, in addition to specifying the inflow internal energy 
by 0, the pressure is assumed to be uniform across the 
lower boundary. Its value ptot is determined to keep the 
mass in the simulation box unaltered (see equations (21) 
to (23) in the cited reference). In downflow regions, the 
velocities and the entropy density S are set according to 

du _ <3u _ dw_ _ ^ £^_q / 14g \ 
dx dx dx ' dx 

(see equations (16) a nd (17) in the same reference). 

In upflow regions, IVogler et al. (<2005l) force the inflow 
to be vertical (cf. (18) therein). In contrast, 



du dv dw 
dx dx dx 







(14b) 



in our implementation. The specific internal energy of an 
inflow is found by ([7]) already discussed further above. 

Open bottom boundaries can lead to a strong increase 
in horizontal momenta since, due to (O, they do not con- 
serve mass and momentum simultaneously. To mitigate 
this effect we usually damp the horizontal momenta with 
the procedure described in Scction l2.3.3l in the layers where 
the bottom boundary condition is applied. 

Table [1] provides a summary of models and boundary 
conditions which we test in Section [ 



2.3. Simulation Domain Size and Initial Conditions 

Numerical simulations of stellar surface convection are 
commonly initialised with one-dimensional models taken 
from stellar evolution calculations. These models often 
do not contain enough data points in the photosphere or 
introduce other deficiencies, which do not allow to apply 
them directly. In general we have to extend the starting 
models both at the top and the bottom. In the following, 
we describe the procedures used for this purpose. 

2.3.1. Extending the Simulation Box in the Photosphere 

To decrease the influence of the top boundary on the 
simulation, it is desirable to move the top boundary away 
from the optical surface towards the upper photosphere 
or even further above. Since the one-dimensional mod- 
els used as an initial condition often do not contain data 
for sufficiently low optical depth and moreover since tur- 
bulent pressure and multidimensional radiative transfer 
cause an elevation of the opt ical surface in comparsion 
with one-dimen s ional models ( Nordlund and Steinl . 1999; 
Rosenthal et al. , Il999h . we developed a procedure to ex- 
tend the simulation domain when setting up the initial 
conditions from such models. 

The basic variables are re-allocated to the new number 
of grid points and are pointwise continued to the (new) 
top of the computational domain by extrapolating the ex- 
isting values. First of all, the pressure p is exponentially 
extrapolated from the values in the initial model. Then, 
the density p is set according to the equation of hydrostatic 
equilibrium 



dp 
dx 



-P9, 



(15) 



where the derivative of p is calculated numerically with 
fourth order centred differences (a negative sign for g com- 
pensates for x increasing with depth). By calling the equa- 
tion of state, T = T(p, p) is calculated. With the new 
values of T and p the corresponding pressure values are 
calculated again and the whole procedure is repeated sev- 
eral times to ensure consistency with the equation of state. 
The velocities are then set according to (15a)) . 

In any case, the initial model for the new layers is only 
of minor importance, since due to the short time scales in 
the upper photosphere, the state variables change rapidly 
away from the initial values. It is indeed more important 
that this procedure is sufficiently robust. The additional 
space allows the model to relax to cooler temperatures and 
reach smaller optical depths. 

2.3.2. Deep Models 

When we tried to calculate deeper solar granulation 
models which extend to a depth of about 4 Mm below the 
surface, we experienced stability problems when the val- 
ues from the one-dimensional starting model were used for 
initialisation of the state variables. This probably stems 
from the fact that the initial model is not numerically adi- 
abatic and hydrostatic when interpolated to the numerical 
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name 


type 


detailed description 


equations 


time scales [h] parameters 


BC 1 
BC 2 
BC 3a 
BC 3b 


closed 
open 
open 
open 


Muthsam et al. f 2010a) 
Voder et al. (2005) 

similar to Frevtag et al. (2012) 


a 

@, © & (|l4b| 

©, ©, CP & d 
©, HDD, GO & CCD 


tkh ~ 550 

T S W 1000 C Pchangc = 1.0 
T S W 100 C Pchangc = 0.1 



Table 1: Summary of bottom boundary conditions in ANTARES. In the fourth column, the number of the equations where the boundary 
conditions are described in this paper are given. 



grid used in the simulation. The deviations from equilib- 
rium and the ensuing instabi li ties g et larger the deeper 

experienced similar 



Robinson et al 



the model is. 
problems. 

Therefore, in the solar case we re-integrate such initial 
models starting from around 2.5 Mm below the surface us- 
ing the classical, explicit 4 th order Runge-Kutta algorithm 
to solve the equations of adiabaticity and hydrostatic equi- 
librium, i.e. 



dp 

to = ~ A9 ' 
dT T 
-rr~ = ~pg— v ad , 
ox p 



(16) 
(17) 



with V a d and p for given p and T taken from the equation 
of state. 

Since we have assumed the turbulent pressure to be neg- 
ligibly small here, the open bottom boundary has to be 
located sufficiently deep below the superadiabatic layer, 
where pturb reaches its maximum ins i de th e convection 
zone (cf., for instance, iTanner et al.1 I2012T) . For very 



deep extrapolations of more than about two pressure scale 
heights, S'inflow can no longer be considered a constant, 
since it will slightly increase even inside the nearly adia- 
batic part of the convection zone, in which case necessary 
changes to S'inflow would have to be estimated from solar 
(resp. stellar) structure models. 

With this procedure a stable initialisation of simulations 
of deep stellar surface convection such as that one of our 
Sun is readily possible. 

2.3.3. Damping Initial Oscillations 

Since initial models set up from one-dimensional stellar 
models numerically are not in perfect hydrostatic equilib- 
rium, every multi-dimensional model starts to oscillate in 
the vertical direction. These oscillations are not the p- 
modc oscillations expected to be seen in the simulation of 
stellar surface convection, but are an artifact introduced 
by the setup of the initial conditions and must be removed. 
Therefore, we adapt ed t he procedure described in d etail in 
iTrampedach (1997) and Trampedach et al. ( 2013 ) in our 
code to damp unwanted vertical oscillations. The damping 
term is an additional term in the velocity equation. The 
damping velocity u mo dc was defined by 



^modc 



(P u ) 
(P) 



(18) 



The (vertical) momentum equation takes the form 

9P U , . "mode , 1fl x 

— - — = usual terms — p , (19) 

ot ^modc 

where £ mo dc is the period of the mode which should be 
damped. To be consistent, we also added a corresponding 
term to the total energy equation such that 



9E fmodc 

-—- = usual terms — pu . 

(st ^modc 



(20) 



which improved the stability of this procedure and re- 
moved its non-conservativeness with respect to energy 
(even though the conservation errors of the original 
method are very small). 

This procedure can easily be extended to remove dis- 
pensable horizontal momenta by replacing the vertical ve- 
locity in the calculation of u mo de by its horizontal counter- 
parts f mo de,ij Or f mo de,2 : 



^modc.y 



(P) ' 



(p) 



(21) 



The time scale i mo dc can be used to control the rate at 
which the horizontal momenta are removed. With values 
of tmodc ~ 1 sound crossing time, we found that removing 
of horizontal momenta, if desired, works very efficiently. 

2.4- Asymmetric Stencils at the Domain Boundaries 
As described in lFerziger and Peric ( 2002 . p. 53), enforc- 



ing conditions which set derivatives to zero as in various 
prescriptions for the in- and outflow at top and bottom 
boundaries like ([5"a|). ©, (|14al) . or (I14b[) should not be 
done by simply setting u to a constant value, but by ex- 
act inversion of the stencils used for the calculation of the 
derivatives. Especially, if higher-orde r methods as, e.g. 



the fourth order method presented in lHappenhofer et al 



(l2013h . are used, unsuitable procedures can lead to patho 
logical behaviour of the velocity field near the boundary. 

In general, boundary stencils cannot keep the order and 
the width of the interior stencils at the same time espe- 
cially when building higher order derivatives. Either the 
stencils are broadened to keep the order of the method con- 
stant, or the width of the stencils is kept constant leading 
to a lower-order method. We decided to conserve the or- 
der of the method and present the corresponding stencils 
in the following. The symmetry of the stencils gets lost in 
any case. 



Boundary conditions for ANTARES must be set on 
three vertical layers. All physical layers are numbered 
from 1 to n x starting at the top of the simulation box, 
such that the boundary layers at the top of the domain 
get the indices 0,-1 and —2, whereas n x + 1, n x + 2 and 
n x + 3 correspond to the boundary layers at the bottom. 

At the top boundary of the computational domain, 
derivatives of a scalar-valued function at the cell cen- 
tres are calculated by 



-250- 



- 360o + 160i — 30 2 



dx 



12Ax 

-30_2 - 100-1 + 



+ 02 



12Aa 



(22a) 
(22b) 



follow ing the approach presented in lHappenhofer et al 
( 20131 ). For |^(xo), the symmetric stencil can be used 
i.e. 



Tx^ 



-l +i 



12Ax 



(22c) 



where 0j = <j)(xi) and mirrored expressions are valid for 
the bottom boundary. For §j(xi) and other derivatives 
at the cell boundary, similar expressions can be calculated 
with the same procedure. 

Expansion in Taylor series shows that the error term e 
is given by 



-i(Ax) 4 0( 5 >(C), (23a) 
5 

e (^(z_0)=^(Ax)V 5) (C), (23b) 

£ (g(*o))=-^(A.)V 5) (C), (23c) 



£ ' Tx {x - 2) 



demonstrating the fourth order of the given stencils. The 
error constant increases the more the stencils are asym- 
metric. 

If we try to enforce §j = by setting 0_2 = 4>-i — 0o = 
0i, the derivatives will not be zero numerically. Instead, 
we reset according to 



(24a) 
(24b) 
(24c) 



Plugging equations ((24]) into (l22l gives 0.0 for the deriva- 
tives numerically, as desired. 

While this procedure immediately applies to bottom 
boundary conditions as specified through (J7J, ©, (|14a[) . 
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- 2 =55 C 
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63 






- 1= 55 C 
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64 






^°=55 C 




h 



and (|14bp . for conditions set by ([6]), ([9]) or (ITUt in con- 
junction with (TlUl) . the vertical velocity u in the innermost 
boundary layer n x + 1 is already fixed by the enforcement 
of mass conservation (| 13d|) . Whereas we can use the proce- 
dure from above for the horizontal velocities, the boundary 
condition for the vertical velocity must be changed to 



u„ x+ i set to enforce mass conservation, (25a) 
279 99 17 

u».+2 =Y97 U "* +1 ~ T97 U "^ + 197""-- 1 ' (25b) 

252 64 9 . . 

Un.z+3 =Y^U„ x + i - J^Un x + jQ = Un x -l. (25c) 

2.5. Calculation of the Energy Fluxes 

The calculation of the energy fluxes follows ICanuto 
(|1997j ) which describes the procedure for the case of an 
ideal gas. We avoid simplifications in the equation of state 
in the following, and set 



(P) 



(26) 



The equivalent holds for v" and w". A positive flux 
means flux towards the top of the box and vice versa. 

2.5.1. Radiative Flux (F ia d) 

In optically thin regions, the radiative transfer equation 
is solved with the short cha racteristics method as described 
in Muthsam et al. ( 2010aF) . In this approach, the equation 



is solved by an angular integration over a discrete set of 
rays. The dependence on the frequency is considered by 
repeating the calculation several times for averages over 
frequency sets, to each of which the wei g ht o^bm is assigned 
( Nordlundl . 119821 IStein and Nordlundl . 120031) . In every 
computational node, the intensity /(bin, ray) is calculated 
for Arrays rays, where A rays depends on the quadrature 
formula chosen for the angular integration , and Abms fre- 
quency bins. For the integration rule from Carlson (1963) 
which is used in the simulations presented in this paper, 
Arays is 24 in three and 12 in two dimensions due to sym- 
metries. Abins = 1 corresponds to the grey approximation, 
and Abins = 4 is used for non-grey simulations. 

The vertical component of the radiative flux is calcu- 
lated by 



:ad 



47T 



bins 



rays 



n x J(bin,ray), (27) 



where n x is the vertical component of the ray direction. 
In optically thick regions, the diffusion approximation is 
valid and F ra d can be calculated simply by 



-Frad = « 



dT 
dx ' 



(28) 



The sign in front of the thermal diffusivity n is positive 
due to the coordinate system chosen in ANTARES. 
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In most of our simulations, the transition to the diffu- 
sion approximation is done at a fixed geometrical depth. 
We choose this depth such that the transition is in any 
case performed in an optically thick region, also in case of 
optically thin downflows. 

2.5.2. Convective Flux (F conv ) 

The vertical component of the convective flux is defined 

by 



F com = (u"p (h-£±£\ 



(29) 



where h = ^-y- is the specific enthalpy. 



2.5.3. Kinetic Flux (F kin ) 

The vertical component of the kinetic flux in three di- 
mensions is defined by 



F kin = ^{pu" (i 



(30) 



When the model is statisti cally stationary, thes e fluxe s 
coincide with those defined in lNordlund and Steinl (|200ll ). 
Furthermore, the viscous flux is defined as the horizon- 
tal avera ge of the viscous st r ess te nsor as in equation 
(16) from iNordlund and Steinl ([20011 ) or equation (7b) in 
Canutol (|l997h . 



3. Results 

In the following, we show results from a bunch of simu- 
lations of solar surface convection with ANTARES. With 
these data, we demonstrate the advantages and disadvan- 
tages of the methods presented in Section [2] 

These results can be generalised to other main sequence 
stars without any fundamental complications. Neverthe- 
less, care must be taken in choosing the boundary condi- 
tion parameters adequately. Since this paper is focussed 
on validating the methods from Section f2J we restrict our- 
selves in the following to the solar case as the best-known 
example of stellar surface convection. 

3.1. Top Boundary Conditions 

There is one free parameter in the top boundary condi- 
tions described in Section \'2. 2.21 Changing the parameter 
Cf controls the speed at which (e)o, the value of the specific 
internal energy in the upper boundary layers, is adapted 
to the values of e in the inner domain. 

In the solar case, typical values of the mean speed of 
sound (v sn d)i and the mean pressure scale height (H)\ 
in the uppermost domain layer are 7km/s and 100 km, 
respectively. Given a time step size per Runge-Kutta in- 
tegration ofjuound 0.06 s, the value of 10~ 3 given for 5 in 
Cheund (l2006r) corresponds to Cf ~ 4.0. 

In two-dimensional numerical experiments we found 
that the value Cf = 4.0 leads to a very slow change in 
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Figure 3: Vortical profiles of at, quantifying the temporal variation 
of the mean temperature (see I31I 0. obtained for different values of 
Cf as defined in J5c tt . 



(e)o- The (temporal) standard deviation of the mean tem- 
perature stratification cr t , 



a t (T) (x) — (((T)(x) - ((T)) t (x)) 2 )t 



(31) 



where (T) is the horizontal and ((T)) t the temporal and 
horizontal average of the temperature T, is plotted in Fig- 
ure [21 The rightmost three nodes in this figure correspond 
to the boundary layers. A value of Cf = 4.0 decreases the 
variation to an undesirably small level at the top, which 
eventually leads to an unrealistic bump in the average T(t) 
relation for layers near the top boundary due to a forced 
slow relaxation, whereas Cf = 0.1 leads to much stronger 
variations in the same region. We suggest to use a value 
of Cf « 0.4 or less. Overall, the effect of Cf on the mean 
temperature stratification throughout the model domain 
is quite small and restricted to the upper layers of the 
simulation. Nevertheless, while the choice of large Cf can 
hinder a fast relaxation to a different statistically station- 
ary temperature structure, equation (|5e[) leads to a time 
scale which is independent of the numerical setup of the 
simulation. On the other hand, very small values of Cf may 
lead to strong variations in the temperature which can de- 
crease the numerical stability of the boundary condition. 

The value of Cf = 0.4 suggested here is a compromise be- 
tween fast relaxation and robustness with respect to strong 
shock fronts. We note here that — also for stability rea- 
sons — a very "stiff" boundary resulting from Cf > 1 is not 
desirable either. Given the numerical justification of this 
choice however, it has to be expected that for stars very 
different from the Sun, for instance, Cepheids or A-type 
stars, Cf will have to be adjusted. 

3.2. Bottom Boundary Conditions 

We compare three 3D solar models with the same sim- 
ulation setup, which differ mainly in the bottom bound- 
ary condition. The simulation boxes extend to 5.2 Mm in 
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the vertical and 9.0 Mm in the horizontal direction. In 
the vertical direction, the resolution is 15.3 km whereas in 
the horizontal directions, we chose a resolution of 32.1 km. 
The term "optical surface" in the following designates the 
position where the mean temperature gradient is steepest. 
This occurs at a geometrical depth of around 0.7 Mm to 
0.8 Mm from the top of the simulation box. That region 
is typically located some 0.1 Mm below the layer where 
TVoss ~ lj but is more easily identified when comparing dif- 
ferent simulations. For the radiative transfer equation, we 
used the grey approximation in the upper 30 % of the box, 
corresponding to a geometrical depth of around 1.6 Mm, 
and the diffusion approximation elsewhere. The time av- 
erage covers about one hour. 

In Model 2, the bo ttom boundary condition BC 2 as 
m IVoeler et all (|2005h was used, with the formula ([7]) for 
calculating the energy of the inflow. Model 3a uses the 
boundary conditions BC 3a and Model 3b uses BC 3b as 
summarised i n Tabled] They are similar to the boundary 
conditions of iFrevtag et al.l ([20121 ). but with different val- 
ues for cp c hange and different formulae for the calculation 
of inflow (Unj is used for Model 3b and © for Model 3a. 
Whereas we use tkh ~ 550 h as time scale in (0 for Mo- 
del 2, we chose r s « 100 h for Model 3b and r s 1000 h 
for Model 3a. Of course, these time scales are not univer- 
sally applicable, but even for solar convection simulations 
have to be adjusted to the depth of the simulation box. 

3.2.1. Energy Fluxes 

The radiative flux at the top of the domain depends 
on the amount of energy flowing through the simulation 
domain. The idea of formula (O therefore is to correct the 
(unknown) entropy of the inflow at the bottom boundary 
according to the radiative flux at the top of the domain 
until the radiative flux has reached the desired value. 

Equation (20) in lVogler etail (|2005l) . i.e. © in this pa- 
per, defines a similar procedure and corrects the specific 
internal energy e at the bottom boundary on the Kelvin- 
Helmholtz time scale tkh ■ Tkh increases rapidly with box 
depth, as shown in Figure [5] For our standard simulations 
of solar surface granulation with a domain reaching two to 
four Mm below the optical surface, tkh ~ 100 h to 600 h. 
In contrast, the time scale for the surface granulation is 
about several minutes. 

In Figure^ the radiative flux at the top of the domain of 
Model 3a and Sinflow calculated with © are plotted. Due 
to the insufficient vertical resolution of the simulation of 
only 15.3 km and the short relaxation time from the one- 
dimensional model, among others, the radiative flux stays 
several per cent too high. 

Figures [5] and [HI show the time-averaged convective and 
kinetic fluxes of Model 2, Model 3a, and Model 3b from 
above, with the fluxes computed by the formulae (|29p and 
([50")l . respectively. 

It is obvious from Figure |4] that the correction at the 
bottom boundary does not influence the radiative flux at 
the top of the domain even after more than two hours 
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Figure 4: Effect of using a variable inflow entropy as in J9} for cal- 
culating Si n fl ow in Model 3a. Radiative flux F^°^ calculated as hori- 
zontal average of F ra[ j in the outermost domain layer and normalised 



of simulation time. During this time, the entropy of the 
inflowing material is decreased by about 0.4% which is 
reflected in a completely wrong overall flux at the bottom 
boundary shown in Figures [5] and El This implies that the 
correction will also not work if applied on even shorter time 
scales than tkh because then, the unfavourable corrections 
due to ([H]) would be even larger. 

We deduce that the coupling of inflowing energy at the 
bottom boundary and the energy radiated at the top works 
on much longer time scales than common in simulations 
of solar and stellar surface granulation. The formulae (JJJ 
and © require a much longer pre-relaxation time of the 
model to be effective, or a model which is already relaxed 
to a suitable three-dimensional stratification. 

From Figure [5] we also deduce that blindly employing 
formula © as used in BC 2 led to a completely wrong 
convective flux at the bottom boundary. The flux is too 
large by a factor of 4. In c ontrast, the energy fluxes from 
IStein and Nordlundl (l2000h do not show any anomaly near 
the bottom boundary. We suppose th at the (fixed in time) 
value for the energy of an inflow in Stein and Nordlundl 
(2000J) was well-chosen, whereas we started with the value 
from the one-dimensional initial model which cannot be 
expected to be in agreement with the value suited for a 
three-dimensional simulation, and which in the end in- 
voked an unreasonably strong correction on iSmflow 

Due to the insufficient vertical resolution of the sim- 
ulation, the radiative flux at the top does not approach 
the desired value no matter how strong the correction at 
the bottom boundary is — at least not on time scales af- 
fordable for a three-dimensional simulation. Since tkh is 
that large and the coupling between the radiative flux at 
the top and ei n fl ow or Sinflow too weak, the convective flux 
is forced to unacceptable values by the bottom boundary 
condition BC 2. Similar is valid for Model 3a where (JSJ) is 
used. There, -Fconv even changes sign. 
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Figure 5: Convective flux of three solar models which differ in the 
bottom boundary condition, normalised by F+ = crT^. The mean 
pressure profile from Model 2 is shown, too (the pressure profiles 
of the models do not differ significantly). The boundary conditions 
affect the lower 2 Mm, or two pressure scale heights, of the domain. 



This weak coupling is also demonstrated by the fact that 
the radiative fluxes of all three models do not differ signifi- 
cantly, even though the differences in the lower part of the 
simulation domain shown in Figures [5] and [5] are exorbi- 
tant. To be more precise, in the lower 2 Mm of the models, 
or two pressure scale heights, the fluxes are influenced in a 
very direct and strong way by the boundary condition de- 
spite their main difference is the determination and change 
of ei n fl ow (resp. Sinflow) as a function of time (Models 3a 
and 3b also differ in the damping rate of pressure fluctua- 
tions at the bottom boundary, cp cria ngc)- Near the surface 
however, the differences between the three models are in- 
significant. 

We note here that the simulations presented in 
Voder et al. have much more shallow domains with 

a vertical extent of typically just 1.4 Mm of which only 
0.8 Mm are located below a mean depth for which the op- 
tical depth is 1. As a result, in their case tkh ~ 1 to 
2 h. For this special case we can expect a much tighter 
coupling of bottom boundary and F T ££ than we observe 
here for our much deeper models because the distance be- 
tween the bottom boundary and the optical surface is only 
about two pressure scale heights in which case a correction 
procedure as given by ([7]) can still be effective. 

We hence prefer BC 3b as used in Model 3b, since these 
boundary conditions keep the energy fluxes in the lower 
part of the domain at the desired level and do not produce 
any serious artifacts even if the initial value of Sinflow is not 
chosen very well. The variation over one time step is small 
for this procedure. If desired, a (practically) constant in- 
put entropy can be obtained by choosing rg tkh- 

Since we do not use any artificial diffusive terms in our 
simulations presented here, the mean viscous flow is very 
small such that we decided not to include it in our plots. 
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Figure 6: Kinetic flux of three solar models which differ in the bottom 
boundary condition, normalised by F+ = crT^g. 



3.2.2. Dependence on Parameters of the Boundary Con- 
ditions 



The parameter cp c han g e in iFrevtag et al.l ([20121 ) defines 
a time scale on which pressure fluctuations at the bottom 
boundary are damped, rs controls the time scale on which 
■Sinflow is changed. We tried two values for both param- 
eters and found that they influence the behaviour of the 
boundary conditions considerably. 

In Figure El the kinetic flux of Model 3a goes to near 
the bottom boundary giving a completely wrong flow pat- 
tern in the lower half of the simulat ion domain. F o r Mo- 
del 3a, the parameter cp cna ngo from IFrevtag et all ( 2012 ) 
was fixed to 1.0 and r s w 1000 h in HTUJ). The only rea- 
sonable fluxes are found in Model 3b, where the formula 
(fTU)) was used for the calculation of Sinflow with rg w 100 h 
and where cp c hange was set to 0.1. It is striking to see how 
large the influence of these parameters can be. Model 2 
uses BC 2. 

In Figure [3 the distribution of the vertical velocity of 
the three models in three vertical layers with fixed depth 
and scaled by the maximum velocity in this layer is de- 
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picted, i.e. the distribution density function of the dimen- 
sionless variable u sca i e d with 
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Figure 7: Scaled velocity distribution 1.3 Mm, 2.8 Mm and 4.3 Mm 
below the optical surface (from top to bottom). The velocity is scaled 
by the maximum velocity in the corresponding layer. For theses pic- 
tures, the velocities in each of the layers were grouped into 96 equal- 
sized bins and scaled by the maximum velocity in this layer. Finally, 
the distribution function was normalised by the number of nodes and 
time steps such that the integral from —1 to 1 is 1. Due to the co- 
ordinate system chosen for ANTARES upflows imply u(x, y, z) < 0. 
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In the layer corresponding to the uppermost image in 
Figure [3 which is situated 1.3 Mm below the optical sur- 
face, the differences between the three models are very 
small. As expected, the distribution has a significant skew- 
ness due to the asymmetry between the faster, narrow 
downflows and the slower upflows. Note the shift of the 
maximum induced by the broad upflows (and the fully 
compressible treatment of the flow). Going deeper and 
therefore closer to the bottom boundary the differences 
between the models become more prominent. In the bot- 
tom panel, the distribution of Model 3a has much narrower 
tails and its skewness is much smaller than in Model 2 or 
Model 3b, whereas the velocity distribution of Model 2 and 
Model 3b is in all three layers very similar (as we expect 
it to be). On the other hand, there are some deviations 
visible for Model 2 while for Model 3b, if scaled to the 
same maximum value, i.e. 



calcd 



[x,y,z) 



u scaied (x, y, z) 
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the distribution is almost invariant as a function of depth 
apart from a gradual broadening observed for layers far- 
ther away from the optical surface. By modifying only the 
correction formula for Si„j ow and two parameters in the 
bottom boundary conditions, we changed the velocity dis- 
tribution of Model 3a in the lower third of the simulation 
domain considerably: the distribution approaches a nar- 
row Gaussian one with a very flat, additional tail and, in 
the end, resulting in a kinetic energy flux of a bout 0. 



Com paring our results with the data from IStein et al 



who performed simulations of the upper part of 
the solar convection zone in boxes of up to 20 Mm depth, 
we conclude that there is no physical reason for the kinetic 
flux to go to 0. On the contrary, it should even slightly 
increase in magnitude when going deeper down into the 
convection zone — even though that trend in the kinetic 
flux in the lower part of their simulation box may also 
be influenced by the boundary conditions. Anyway, the 
kinetic fluxes of Models 2 and 3b are much more realistic 

than the one of Model 3a. 

Fur thermore, we conclude from Figure 3 in IStein et~al 



( 2009f ). that the assymetry between up- and downflows 
persists when going deeper into the convection zone, while 
simultaneously the dominant horizontal scale of the con- 
vection increases. This means that the velocity distribu- 
tion function in each vertical layer should have the same 
shape when scaled to the maximum velocity value in this 
layer. 

Returning to the differences between Model 2 and Mo- 
del 3b, in the bottom panel of Figure [JJ the peak of the 
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scaled velocity distribution of Model 2 is close to a veloc- 
ity of whereas for Model 3b, the skewness is not changed 
significantly. The same can also be observed when com- 
paring more shallow models such as those summarised in 
Tabled What we obse r ve for Model 3b is hence in agree- 
ment with Stein et al. (2009), which is not the case for 
Models 2 and 3a as their behaviour as a function of depth 
may be modified by just changing the depth of the simula- 
tion box. This makes no physical sense, as we require the 
open bottom boundary condition to influence the flow as 
little as possible. Altogether, this indicates that Model 3b, 
and therefore BC 3b is superior. 

Since both values tested for rs are rather large, we as- 
sume that they do not influence the results as does the 
change in the value of cp c hange on the one hand and the 
change from (O to (fTU| on the other. Individual testing 
of the parameters, however, would have to be done also 
in three dimensions and with sufficient numerical resolu- 
tion, which has not been affordable in terms of computa- 
tion time. Anyway, the combination of rg ~ 100 h and 
cpchange = 0.1 yields sensible results and we choose them 
as new default values for these parameters. 

From the comparison of the fluxes we deduce that the 
domain of influence of the bottom boundary condition in- 
cludes the lower third of the simulation domain, in this 
case about 2 Mm or 1 to 2 p r essure scale heights, coinciding 
with the findings inlKupka| d2009bl p. 95) and suggested in 
Kupka and Robinson! (|2007l ) for the case of closed bound- 
ary conditions. The fluxes in the upper two thirds of the 
domain are nearly identical for all three models, demon- 
strating that at least the surface layers and the superadi- 
abatic region are not strongly influenced by the boundary 
conditions, if the safety margin is large enough. 

3.3. Dependence on Initial Model 

In the following, we compare three solar 3D models 
which differ only in their initial conditions, three differ- 
ent ID solar structure models. The first one, called mod- 
els in the following, is based on a semi-e mpirical model - 
It com bines an MLT model of convection (jBohm Vitensd . 
1958) with aMLT = 1-99 with a semi-empirical model for 
the photosphere. It is our standard initial model and 
was used for all other simulations of solar surface con- 
vection in this paper. The second one is a purely the- 
oretical model. Convection is modelled by means of an 
MLT model with qmlt = 1-766. F inally, the thi r d one 



uses the CGM model of convection (jCanuto et all [1996) 
with «cgm = 0.69. Both theoretical models are calcu- 
lated with the CESAM code and use t he grey Eddington 
approximation for the phot osphere. See Morel! ( 1997 ) and 
Morel and Lebretonl (ElO 8J) for a de scription of the CE- 
SAM code and ISamadi et al.l ( 20061 ) for the implementa- 
tion of the convection model. 

Since the photospheric layers contained in the stellar 
structure of the models computed with the CESAM code 
are much more shallow than modelS, we extended the box 
of all three models such that the computational domain 
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Figure 8: Entropy of the three initial models (modelS, CGM, MLT) 
and of the relaxed numerical simulation as well as tkh as a func- 
tion of depth. The data for the initial models are obtained from an 
average over the first 0.6 s of each three-dimensional model. Since 
the simulations based on these initial conditions all feature similar 
relaxed profiles, only one representative simulation is shown in its 
relaxed state, too. 



has about the same size, i.e. 6 Mm in the horizontal and 
4 Mm in the vertical direction and such that the optical 
surface is at about 1 Mm geometrical depth. All models 
treat radiative transfer with the opacity binning method 
with 4 bins and were relaxed for 80 to 90 min with low spa- 
tial resolution. Then, the following 30 min were used for 
analysis. During this time the resolution was about 11 km 
in the vertical and 35 km in the horizontal directions. All 
models use the open boundary conditions BC 3b at the 
bottom. 

As depicted in Figure [51 the initial entropy profiles dif- 
fer mainly in the layers from 1 Mm below the surface up 
to the surface itself. At the bottom boundary, the dif- 
ferences in entropy are much smaller. Immediately below 
the superadiabatic peak, the three profiles from the re- 
laxed simulations, however, are very similar to the CGM 
initial model. As they do not differ considerably from each 
other, only one profile is shown in Figure [8j Since tkh is 
about 1 to 2 h in this region, the pre-relaxation time was 
long enough to allow the simulations to adjust to a com- 
mon stratification. Near the surface (at about 150 km), 
the simulation results differ from all three one-dimensional 
models. On the other hand, when we look at Figure HI we 
see that V — V a d from the simulation is much closer to the 
MLT and the modelS profiles, a result which may be un- 
expected when looking at the seemingly small differences 
in the entropy profile near its minimum (Figure We 
can understand this result by taking into account the el- 
evation of the optical su rface in the fully developed three 



dimensional simulation ([Nordlund and Steinl . 11999! ) . This 
process moves the region of nearly constant entropy up- 
wards while the averaged superadiabatic gradient remains 
less steep in the very narrow region around its maximum. 
In all three models, once relaxed, the entropy profile 
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Figure 9: V — V a< j of the three initial models and of the numerical 
simulation as discussed in Fig. [8] The data for the initial models are 
obtained from an average over the first 0.6 s of the three-dimensional 
model. The integration rule from [Carlson! dl963T) was used for the 
angular integration in the radiative transfer solver. The logarithmic 
temperature gradient V is computed with respect to the total (i.e. 
including turbulent) pressure. 



as well as the radiative flux at the top do not differ sig- 
nificantly from each other. Thus, differences in the su- 
peradiabatic layer exhibited by the initial models are not 
important, if the simulations are run long enough in terms 
of tkh for that region. From Figure [5] we deduce that 
tkh is about 1 to 2 h in the region where the entropy pro- 
files of the initial models differ, which corresponds to the 
simulation time used for pre-relaxation. 

3.4- Resolution Dependence 

To reach the correct radiative flux at the top the res- 
olution of the simulation must be high enough and the 
initial model must be well-suited. From numerical exper- 
iments we deduce that a vertical resolution of at least 10 
to 1 1 k m is necessary for this purpose, if the integration 
rule of ICarlsonl (|l963h is used for the angular integration. 
These numbers are slightly different depending on the type 
of binning used for the radiative transport, and depend 
also on the quadrature formula used to perform angular 
integration. 

The resolution dependence is demonstrated in Table [2j 
where a bunch of grey and non-grey models are compared 
with respect to the radiative flux at the top of the do- 
main. All of these models differ from the ones mentioned 
above in terms of simulation domain size, duration of the 
simulation and grid resolution. Whereas the grid resolu- 
tion is specified directly in Table [2 the simulation domain 
is always between 3.8 Mm to 4.1 Mm in the vertical and 
6.0 Mm in the horizontal direction. All models are three- 
dimensional and the radiative transfer eq uation is solved 
along 24 rays using the quadrature rule of Carlsonl (|l963h 
for the angular integration, but the number of bins differs. 
The grey approximation corresponds to the case of one 
bin, whereas the non-grey models were calculated with 
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Table 2: Effective temperature calculated from F °^ = crT^j in de- 
pendence of binning method, numerical resolution and boundary con- 
ditions. The observed effective temperature of the Sun is 5777.6 K. 
Boundary conditions are explained in Tabled 



four bins. The bottom boundary conditions are specified 
in Table [2] Finally, the model must reach far enough into 
the photosphere. We found that about 800 km above the 
optical surface are sufficient. 

F °¥ primarily depends on the vertical resolution. In 
the grey case, F*°% is slightly lower than in the non-grey 
case where we used 4 bins. We emphasize that the rate 
of change of with numerical resolu tion depends o n 
the angular quadrature rule. The rule bv lCarlsonl (jl963ti . 
which does not include a vertical ray, requires a very high 
resolution to yield the correct radiative flux (only the mo- 
del with 6.5 km vertical resolution yields a T e ff very close 
to the solar one). 

Conversely, the dependence of F*°^ on the lower bound- 
ary conditions is weak. Table [5] does not show any system- 
atic difference with respect to the bottom boundary con- 
ditions. Changes occur only on much longer time scales. 

We note that the horizontal resolution of all simulations 
presented here is sufficient to get the correct effective tem- 
perature and radiat i ve flux at the top of the domain (cf. 
Asplund et ai1l200oHRobinson et al.ll2003l ). 



3.5. Comparison with Closed Boundary Conditions 

The main effect of an open bottom boundary is that 
it permits a free vertical outflow. In the case of closed 
boundaries, the fluid cannot penetrate the bottom bound- 
ary and energy is transported into the domain only by an 
artificial radiative source term. 

In Figure fTUl the energy fluxes of two similar models, one 
with closed boundaries BC 1 and one with open boundaries 
BC 3b, are compared. The time average extends over two 
hours for the model with closed boundaries, and over one 
hour for the model with open ones. Since the resolution 
of the model with open boundary conditions was higher, 
the data was interpolated to the coarser grid of the model 
with closed boundaries. The resolution changes the radia- 
tive flux near the surface, similarly to what was shown in 
Table [2J thus explaining the step in F conv near a height of 
km. In the following, we only discuss the differences in 
the kinetic and convective fluxes. 

Compared to the model with open boundary conditions 
the convective and the kinetic flux stay smaller in magni- 
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Figure 10: Difference in the energy fluxes of a model with closed and 
a model with open boundary conditions normalised by F+ = crT* ff . 
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Figure 11: Temporal evolution of the normalised total kinetic energy 
for a two- and a three-dimensional simulation. 



tude — keeping in mind their opposite sign. The simula- 
tion with closed boundary conditions is not as well relaxed 
as its counterpart with open boundary conditions, which 
is at least partially due to the insufficient resolution. 

Around 1 Mm above the bottom boundary, the kinetic 
flux from the model with closed boundaries reaches a satu- 
rated state. The differences to the model with open bound- 
aries is very small. On the contrary, the effect of the closed 
boundary on the convective flux is much larger. In about 
half of the simulation box, the convective flux differs con- 
siderably from the one with open boundaries (besides from 
the general offset stemming from the insufficient vertical 
resolution). The difference between the convective fluxes 
exhibits a slope for about 2 Mm starting from the bot- 
tom of the simulation box. The closed boundary condi- 
tion hence also influences about 2 pressure scale heights, 
as is shown in Figure [5] for the case of open boundary con- 
ditions (Model 3a). However, in addition, at the bottom 
boundary itself, the stratification of the model with closed 
boundaries is stable and energy is trans ported into the 



domai n by an artificial source term <|4b[) (jMuthsam et al 



2010aL cf. also Robinson et al. 2003 ) which explains the 
behaviour of the fluxes near the bottom. 

3.6. Two-Dimensional Models 

As demonstrated in Figures 2] and [5J thermal relaxation 
of a simulation takes place on very long time scales, compa- 
rable in magnitude u p to the local Kelvin-Helmholtz time 
scale (jKupka . 2009b) unless the mean thermal structure is 
known in advance (we return to this issue below). Thus, if 
a good initial model is not available, a three-dimensional 
simulation may not be run long enough with reasonable 
effort such that it can relax properly. 

Since the computational costs for a two-dimensional 
model are much lower, one is tempted to relax the mo- 
del in two dimensions first and switch to three dimensions 
as soon as the correct amount of energy is contained in 
the simulation domain. Unfortunately, we have found that 



two- and three-dimensional models of solar surface con- 
vection behave quite differently during relaxation. 

In Figure [TTJ the time evolution of the total kinetic en- 
ergy -Ekin scaled by the sum of kinetic and thermal energy 
of a two- and a three-dimensional simulation both start- 
ing from the same one-dimensional model is shown. In 
theory, the kinetic energy should saturate as soon as the 
convection is fully developed, and the simulation should 
reach a quasi-stationary state. Obviously, this happens 
quite rapidly in the three-dimensional case, whereas the 
kinetic energy of the two-dimensional model starts to os- 
cillate with a very high amplitude and without showing 
any signs of saturation. 

From Figure [TT] we deduce that a two-dimensional simu- 
lation intended as a precursor of a three-dimensional sim- 
ulation should not be conducted for a long time and can 
at most be used to pre-calculate the few first sound cross- 
ing times of the model. Even then, it can be hard to 
get rid of the artifacts of the two-dimensionality, espe- 
cially if the model reaches deep into the convection zone. 
When we tried to convert a two-dimensional model which 
was run for a long time (around 60 sound crossing times) 
to a three-dimensional one, we encountered serious prob- 
lems destroying the symmetry in the third dimension in 
the deeper layers, and the relaxation to a truly three- 
dimensional state was not shorter than when the model 
were started directly from the one-dimensional initial mo- 
del. On the contrary, the model inherits a lot of the de- 
ficiencies of the two-dimensional model such as statistical 
non-stationarity and an overall higher vorticity. Relaxing 
the model to a true three-dimensional state took at least 
as long as when starting from a one-dimensional model. 
This was independent of the type and strength of pertur- 
b ation we used t o break the two-dimensional symmetry. 
(l2000h 



Asplund et al 



pointed out that a two- 
dimensional simulation of solar granulation leads to dif- 
ferent stratifications than a three-dimensional simulation 
and that the results from two-dimensional simulations are 
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qualitatively wrong concerning line profiles and element 
abundances. Furthermore, they emphasized that two- and 
three-dimensional simulations need different values for the 
entropy of the inflowing material (about 2% difference). 
As a consequence, it appeared advisable to us to avoid 
using two-dimensional simulations for testing the bottom 
boundary conditions since the resulting entropy profiles 
will differ from their three-dimensional counterparts. 

Nevertheless, pre-computing a model in two dimensions 
can be useful if the two-dimensional phase does not take 
too long, and if the model is shallow such that the time 
scale of motion at the bottom boundary is still quite short. 
For instance for the solar case, the first hour after starting 
from the one-dimensional initial model may be calculated 
in two dimensions, if the model is not deeper than 4 Mm. 

If the stratification of the one-dimensional model is far 
away from that one of the multi-dimensional model, it 
might take a very long time for the model to relax to the 
new quasi-stationary state. In this case, the horizontal 
and temporal average of a long two-dimensional simula- 
tion could be used as a one-dimensional starting model 
for the three-dimensional simulation. In this way, the 
problem of breaking the symmetries and high vorticity is 
avoided, provided one accepts to discard with any a priori 
information about the velocity field. 

Figure [11] furthermore demonstrates that the thermal 
stratification of a model is quite robust against kinetic mo- 
tions since the kinetic energy reaches only about 0.1% of 
the thermal energy. This ratio is comparable to the ra- 
tio of T S ound to tkh and thus essentially the time scale of 
convective transport relative to the time scale of simula- 
tion relaxation from our "arbitrary" initial conditions (see 
Sect. 

Confirming the results from Section 13.31 convective mo- 
tions would need a long time to change the thermal 
stratification of the entire model, especially since quasi- 
adiabaticity keeps the sizes of the fluxes of kinetic energy 
and enthalpy within the same order of magnitude. A sig- 
nificantly different stratification could of course lead to a 
higher total flux which in turn could reduce the relaxation 
time somewhat. Still, a good starting model is an essential 
prerequisite for cost efficient numerical simulations. 

In Figure [TO. another difference between two- and three- 
dimensional simulations is shown. There, the relative 
share of the horizontally averaged turbulent pressure in 
the innermost boundary layer of the open top boundary 
condition is plotted. The simulations are the same as in 
Figure [TO In both cases, Cf = 4.0, but we tested several 
different values and found no noticeable variation. The 
difference in the relative contribution of turbulent to to- 
tal pressure is significant: in the three-dimensional case, 
it fluctuates around 10% whereas in two dimensions, the 
turbulent pressure sometimes even exceeds the sum of the 
gas and the radiative pressure. In the mean, its share is 
around 30% demonstrating once more the unstable and 
non-stationary nature of the two-dimensional simulation 
for the upper layers of the (solar) photosphere. 
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Figure 12: Temporal evolution of the ratio Pturb/(j> + Pturb) m the 
innermost boundary layer for the same simulations. 



In both cases, the share of the turbulent pressure in the 
boundary layers is not negligible. Using (|5fj) is therefore 
definitely an improvement compared to using the hydro- 
static equation without turbulent pressure. Nevertheless, 
Figure [TO. demonstrates that any static approximation of 
these layers is rather unrealistic. 

(bottom) in 



In comparison with Figure 1.3 (bottom) in ISteffen 
( 2000h . the difference between two and three dimensions is 
of a similar magnitude, even though the layer depicted in 
Figure II 2 | is e ven h igher in the photosphere than the data 
from lSteffenl (1200(1 



4. Discussion 

The physical realism of simulations of stellar surface 
convection depends strongly on the boundary conditions 
imposed at the top and the bottom of the simulation do- 
main. In a large part of the simulation domain, the flow 
is influenced in a direct way by the boundary conditions. 
Therefore, the simulation box must be wide enough in all 
directions, and the boundary should be open to allow free 
in- and outflow. The influence of the lower boundary con- 
ditions is particularly strong on the lowermost two pressure 
scale heights. With closed boundary conditions, the flow 
pattern is disturbed such that t he results are physically not 
relevant in thi s region (cf. also lKupka and Robinson1l2007l : 
iKupkal l2009bT ). Open boundary conditions are clearly 
preferable as they can increase the realism of the simula- 
tion throughout the entire domain considerably, but even 
with open boundary conditions we recommend to ensure 
a distance of at least two pressure scale heights between 
the lower boundary and the superadiabatic layer. 

At the top of the domain, time scales are very short 
such that the combined influence of the top boundary con- 
ditions and the initial stratification of the layers in that 
region on the stratification of the model throughout most 
of the simulation box is rather small provided that it is 
located sufficiently far above the optical surface (cf. Fig. 1 
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in Kupkal 2009a|) . The boundary conditions © assume 
an isoenergetic (isothermal) and hydrostatic stratification, 
including the effect of turbulent pressure. Since the veloc- 
ities are set constant in the top boundary layers, waves are 
not reflected, but transmitted, which increases the long- 
term stability of the model. The formulation of the top 
boundary conditions allows variations in the temperature 
away from the initial profile. In our implementation, the 
parameter Cf controls the stiffness of the open boundary 
conditions in terms of temporal variations of the temper- 
ature, as shown in Figure [3] The optimal value has to be 
found by numerical experiments and may vary depending 
on the specific application. 

The most difficult part in designing open boundary con- 
ditions at the bottom of the simulation domain is to specify 
the amount of energy of an inflow. This value controls the 
overall stratification of a surface convection model. But 
the coupling with the energy radiated at the top of the 
domain is rather weak and for deep, nearly adiabatic con- 
vection zones it operates on time scales which are not fea- 
sible to be calculated with multi-dimensional simulations. 

In deep stellar convection zones the only mechanism 
which can change the stratification of a model is convec- 
tion. But to change the energy content of the lowest 1 Mm 
of one of our solar models by only 1%, convection will need 
about 1 h, while this time scale increases rapidly with box 
depth. The Kelvin-Hclmholtz time scale tkh may not 
only be interpreted as the time it takes until the whole 
gravitational energy contained in the simulation box is ra- 
diated, but also as the time it takes convection to rear- 
range the complete stratification of a (in our case nearly 
adiabatic) model. This results from the fact that (grav- 
itational) potential energy completely dominates the to- 
tal energy balance. For radiative (and conductive) zones 
a derivation of why tkh can be used as an approxima- 
tion to the thermal adjustment time, which is the actual 
quantitiy of interest here, can b e found in Chapter 5.3 of 
Kippenhahn and Weieertl ([1994I ). where also limitations of 
this approximation are pointed out. For the case of con- 
vection zones the approximation tkh ~ ^adj is derived in 
Chapter 6.4 of the same reference as part of a general sta- 
bility analysis for stellar material at a certain (arbitrary) 
depth. As illustrated with Figure [51 differences in the en- 
tropy of the initial model within the superadiabatic region 
disappear within the local tkh whereas the stratification of 
simulations over shorter relaxation times is determined by 
the initial model. The success of one-dimensional models 
in predicting the nearly adiabatic stratification for most 
but about the upper 1 Mm of the solar convection zone 
is essential for the efficiency of hydrodynamical simula- 
tions and for the short relaxation times of 1 to 2 hours of 
numerical simulation of solar granulation reported in the 
literature. 

For simulations of solar convection, any correction mech- 
anism which correlates the bottom boundary directly with 
the top boundary will be inefficient, if the simulation box 
spans more than 1 or 2 Mm of depth below the optical sur- 



face. On the other hand, in shallow simulation boxes the 
boundary conditions will directly influence the superadia- 
batic layer. The simulation box must be large to keep the 
influence of the bottom boundary small. This, however, 
prohibits a direct correction mechanism of the entropy at 
the bottom boundary based on the radiative flux on the 
top. Consequently, the input entropy (resp. internal en- 
ergy) is much better controlled by a local criterion like 
(TTU)) . We prefer to impose the boundary conditions on the 
entropy since its value is nearly constant with depth, such 
that the boundary conditions can be applied at a range of 
geometrical depths using the same parameter values. Cri- 
teria like (J7J or © which couple the upper boundary to the 
lower boundary can lead to wrong energy fluxes as shown 
in Figures [5] and [HI since the correction at the bottom does 
not affect the upper boundary in the time covered by the 
simulation. In terms of Table [TJ our preferred boundary 
conditions are BC 3b. These boundary conditions allow 
free in- and outflow with an adiabatic inflow, while fixing 
the mean mass flux over the bottom boundary to 0. We 
refrain from using BC 2 and BC 3a for future simulations. 

In the formulation of open boundary conditions, several 
parameters must be set which influence the behaviour of 
the flow near the physical boundaries considerably. But, 
since the computational costs for each single simulation of 
surface convection are high, systematic tests of the specific 
effect of each parameter are expensive. Nevertheless, we 
found a reliable choice for all parameters at least for the 
solar case, which should be widely applicable to stars of 
the central and lower main-sequence. 

The computational costs cannot be reduced by perform- 
ing simulations only in two dimens ions, since the result s 
differ considerably. As stated also in lAsplund et al.l ( 2000h . 
two-dimensional simulations lead to different stratifica- 
tions because of the different nature of the kinetics of the 
model. Furthermore, longer simulation times are needed 
for statistical significance, since the area coverage of two- 
dimensional simulations is smaller — which dilutes the 
gain in computation time again. We have found that in 
addition two-dimensional simulations of solar surface con- 
vection do not reach an energetic equilibrium on the same 
time scale as three-dimensional ones and the upper pho- 
tosphere is influenced by much stronger fluctuations of 
turbulent pressure. As the symmetries in the flow-field 
are also difficult to break by applying a perturbation to 
initialise a three-dimensional simulation, the use of two- 
dimensional models as precursors of three-dimensional 
ones should hence be limited to only five to ten sound 
crossing times. Alternatively, horizontal averages may be 
used to construct one-dimensional initial models, if the 
latter cannot be constructed in a reliable way otherwise. 

The numerical resolution influences directly the radia- 
tive flux at the top of the domain. As we show by nu- 
merical experiments in Table [21 the accuracy of the ra- 
diative transfer solver can be improved by increasing the 
numerical resolution on the hydrodynamical grid. If the 
resolution is sufficiently high, the radiative flux and there- 
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fore the effective temperature is in agreement with the 
observed valu es. When the angular quadrature rule by 
Carlson! ([19631 ) is used, a vertical grid spacing of at least 10 
to 11 km is necessary for the solar case, otherwise the effec- 
tive temperature of the model will be too high. Therefore, 
simulations with lower numerical resolution will behave 
differently, and a realistic simulation of surface granulation 
must de done in sufficiently high resolution. If the aspect 
ratio of each computational cell is very different from 1, 
the magnitude of the numerical viscosity varies in depen- 
dence of the coordinate direction, and surfaces which are 
inclined with respect to the coordinate system are not well 
resolved. Therefore, increasing the vertical resolution must 
be accompanied by an increase in horizontal resolution as 
well, which makes the simulation much more expensive. 

The required resolution depends on the angular in- 
tegration rule use d in the radiative transfer solver. 
Tanner et all (|2012h discussed that the local radiative flux 
(in a fixed grid point) is highly variable in time and in the 
inclination angle. Therefore, choosing a suitable integra- 
tion rule for the angular dependence of the intensity func- 
tion is crucial to accurately determine the radiative flux, 
and different formulae will give different results. Only in 
the limit of more and more ray directions these results will 
converge. 

In contras t to the ang ular integration rul e from 
Tanner et ail (|2012h . but also lStein and Nordlundl (l2003h . 
the integration rule used in our simulations (|Carlsonl .ll963) 
does not have a vertical ray. But at least in the centre 
of a granule, the radiative flux is dominated by the ver- 
tical direction. If only few rays are used in the angular 
integration, as it is common in stellar surface convection 
simulations, a vertical ray should be included in the inte- 
gration scheme. We will show in a forthcoming paper, how 
choosing an angular integration rule changes the radiative 
flux of a numerical simulation of stellar surface convection. 
After all, the effective temperature (as well as its supera- 
diabicity) of a model depends primarily on the radiative 
transfer at the surface and only weakly on changes of the 
entropy at a bottom boundary far away from the surface. 

5. Conclusions 

The design of open boundary conditions has to be done 
with great care, especially for high-order methods, and 
their influence on the simulation results can be large. In 
this paper, we investigate the shortcomings of several ap- 
proaches and discuss how to implement boundary condi- 
tions for a high-order method which allow realistic simu- 
lations of stellar surface granulation. 

In a large part of the simulation domain, the properties 
of the flow are strongly influenced by the boundary condi- 
tions for which closed boundary conditions give unphysical 
flow patterns. Hence, open boundary conditions are desir- 
able, except for simulations made for specific applications 
where the flow structure is of limited interest, such as the 
determination of V — V a( j in the superadiabatic layer. 



At the top of the domain, open boundary conditions can 
transmit waves originating at the stellar surface whereas 
closed boundaries reflect them. At the bottom boundary, 
the value of the inflowing entropy must be specified for 
open boundaries and it controls the overall stratification 
in the simulation box. Some of the proposed approaches, 
which couple the input flux at the bottom bound ary to 
Fmd at the top of the domain ( Vogler et all 120051) , work 
only for shallow boxes where the Kelvin-Helmholtz time 
scale is small. Instead, F ra a at the surface of a model is sen- 
sitive to the numerical resolu tion and, as was discussed re- 
cently in lTanner et al. (2012), to the direction and number 
of rays used in the radiative transfer solver. Consequently, 
the entropy of the inflow at the bottom boundary should 
be controlled by a local criterion. In a forthcoming paper, 
we will demonstrate the difference in radiative flux result- 
ing from using different formulae for angular integration 
in the radiative transfer solver. 

Due to the large influence of the boundary conditions, it 
is important to choose the simulation domain sufficiently 
large in all directions. We show how to extend the box if 
the one-dimensional initial models do not provide enough 
data, whence the simulations can be initialised in a nu- 
merically robust way. Great care also has to be taken in 
how to specify either the density, pressure, or momentum 
of the inflow at the bottom boundary, since this can mod- 
ify the fluxes and even the flow topology in the convection 
zone to the same amount as it is done by closed boundary 
conditions. Open bottom boundary conditions can be ex- 
pected to excel over closed ones only, if they are designed 
along these principles and combined with a suitable, local 
criterion for the inflow entropy. We have discussed one 
such successful setup called BC 3b in Section 13.21 

Numerical simulations of stellar surface convection are 
expensive in terms of computation time, since they are 
intrinsically three-dimensional and need high resolution. 
In two dimensions, kinetics are qualitatively different, as 
shown in Section r3.6lwhich extends earlier r esults on strat- 
ification discussed in Asplund et al.l ( 2000l) . Finally, the 
simulations have to be run for a long time to relax to an 
intrinsically three-dimensional state. Since they must not 
be influenced any more by the initial model, this limits 
the usefulness of two-dimensional surface convection sim- 
ulations of the Sun and solar-like stars as precursors of 
three-dimensional ones to short initialisation runs and to 
the computation of mean structure models useful when no 
reliable one-dimensional models are available. 

The high computational requirements of each single sim- 
ulation of stellar surface convection prohibits systematic 
testing of each parameter in the formulation of the open 
boundary conditions, even though the behaviour of the 
latter depends sensitively on them. A larger range of pa- 
rameters and configurations of boundary conditions could 
also be investigated by comparing the results from dif- 
ferent simulations codes. Nevertheless, we have found a 
reliable choice for both top and bottom boundaries while 
we have demonstrated that a poor choice can have a severe 
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impact on a large part of the simulation domain. The best 
performing setup presented in this paper has become the 
default for numerical simulations of stellar surface convec- 
tion with ANTARES. 
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